I'm trying to write a matlab function that creates random, smooth trajectories in a square of finite side length. Here is my current attempt at such a procedure:

function [] = drawroutes( SideLength, v, t)
%DRAWROUTES Summary of this function goes here
%   Detailed explanation goes here

%Some parameters intended to help help keep the particles in the box
RandAccel=.01;
ConservAccel=0;
speedlimit=.1;
G=10^(-8);
%
%Initialize Matrices
Ax=zeros(v,10*t);
Ay=Ax;
vx=Ax;
vy=Ax;
x=Ax;
y=Ax;
sx=zeros(v,1);
sy=zeros(v,1);
%
%Define initial position in square
x(:,1)=SideLength*.15*ones(v,1)+(SideLength*.7)*rand(v,1);
y(:,1)=SideLength*.15*ones(v,1)+(SideLength*.7)*rand(v,1);
%
for i=2:10*t
    %Measure minimum particle distance component wise from boundary
    %for each vehicle
    BorderGravX=[abs(SideLength*ones(v,1)-x(:,i-1)),abs(x(:,i-1))]';
    BorderGravY=[abs(SideLength*ones(v,1)-y(:,i-1)),abs(y(:,i-1))]';
    rx=min(BorderGravX)';
    ry=min(BorderGravY)';
    %
    %Set the sign of the repulsive force
    for k=1:v
    if x(k,i)<.5*SideLength
        sx(k)=1;
    else
        sx(k)=-1;
    end
    if y(k,i)<.5*SideLength
        sy(k)=1;
    else
        sy(k)=-1;
    end
end
%
%Calculate Acceleration w/ random "nudge" and repulive force
Ax(:,i)=ConservAccel*Ax(:,i-1)+RandAccel*(rand(v,1)-.5*ones(v,1))+sx*G./rx.^2;
Ay(:,i)=ConservAccel*Ay(:,i-1)+RandAccel*(rand(v,1)-.5*ones(v,1))+sy*G./ry.^2;
%
%Ad hoc method of trying to slow down particles from jumping outside of
%feasible region
for h=1:v
    if abs(vx(h,i-1)+Ax(h,i))<speedlimit
       vx(h,i)=vx(h,i-1)+Ax(h,i);
    elseif (vx(h,i-1)+Ax(h,i))<-speedlimit
       vx(h,i)=-speedlimit; 
    else
        vx(h,i)=speedlimit;
    end
end
for h=1:v
    if abs(vy(h,i-1)+Ay(h,i))<speedlimit
       vy(h,i)=vy(h,i-1)+Ay(h,i);
    elseif (vy(h,i-1)+Ay(h,i))<-speedlimit
       vy(h,i)=-speedlimit; 
    else
        vy(h,i)=speedlimit;
    end
end
%
%Update position
x(:,i)=x(:,i-1)+(vx(:,i-1)+vx(:,i))/2;
y(:,i)=y(:,i-1)+(vy(:,i-1)+vy(:,1))/2;
%
end
%Plot position
clf;
hold on;
axis([-100,SideLength+100,-100,SideLength+100]);
cc=hsv(v);
for j=1:v
plot(x(j,1),y(j,1),'ko')
plot(x(j,:),y(j,:),'color',cc(j,:))
end
hold off;
%
end

My original plan was to place particles within a square, and move them around by allowing their acceleration in the x and y direction to be governed by a uniformly distributed random variable. To keep the particles within the square, I tried to create a repulsive force that would push the particles away from the boundaries of the square. In practice, the particles tend to leave the desired "feasible" region after a relatively small number of time steps (say, 1000)."

I'd love to hear your suggestions on either modifying my existing code or considering the problem from another perspective.

When reading the code, please don't feel the need to get hung up on any of the ad hoc parameters at the very beginning of the script. They seem to help, but I don't believe any beside the "G" constant should truly be necessary to make this system work.

Here is an example of the current output:

Many of the vehicles have found their way outside of the desired square region, [0,400] X [0,400].